NEUTRON TRANSPORT FROM A POINT SOURCE IN A 
SLAB; A COMPARISON BETWEEN DIFFUSION THEORY 
AND TRANSPORT THEORY 



by 



John Henry Shimerda 




4 «) 



1 



I 



.. I; 

M 

* ’I ^ 



I 



1 




• 1 . 



ii H:! 



-i 



United States 
Naval Postgraduate School 




THESIS 

NEUTRON TRANSPORT FROM A POINT SOURCE 
IN A SLAB; A COMPARISON BETWEEN 
DIFFUSION THEORY AND TRANSPORT THEORY 

by 

John Henry Shimerda 
September 1970 



This document has been approved for public 
release and sale; its distribution is unlimited* 



Neutron Transport from a Point Source 
in a Slab; A Comparison Between 
Diffusion Theory and Transport Theory 



by 



John Henry Shimerda 
Major, UnitecT States Army 
B.S., United States Military Academy, 1958 



Submitted in partial fulfillment of the 
requirements for the degree of 



MASTER OF SCIENCE IN PHYSICS 



from the 

NAVAL POSTGRADUATE SCHOOL 
September 1970 



:ary 

L POSTGRADUATE SCHOOD 
EEEY, CALIF. 93940 



ABSTRACT 

The one-speed, steady state, diffusion equation was solved for a point 
source and for a normal pencil beam in a homogeneous, isotropically 
scattering slab. Numerical results obtained using diffusion theory 
were compared to available transport theory results for two slab thick- 
nesses. 

This comparison demonstrated that the diffusion theory approximation 
to the transport equation will yield accurate results except within 
about one half mean free path of a boundary and except within about 
three mean free paths of a source. The best agreement between diffusion 
theory and transport theory is obtained if the first radial buckling 
constant in the diffusion solution is chosen equal to the radial buckling 
computed using transport theory. 
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r. INTRODUCTION 



After being introduced into a diffusing medium, neutrons suffer 
numerous collisions with the atomic nuclei of the medium. The collision 
process is quite complicated, since either scattering, fission, or 
absorption can occur. As a result of streaming and repeated collisions, 
neutrons are constantly changing their location in the medium as well 
as their velocity (speed and direction). Indeed, some neutrons disappear 
(absorption) while others reappear (fission) while many simply change 
their speed and direction (scattering). The motion of any one neutron 
would appear to be quite random. 

The statistical concept of a distribution function is convenient 
for describing this complicated motion. In essence, we consider "typical" 
neutrons and try to find the neutron density throughout a medium by using 
the principles of transport theory. Although these principles are 
straightforward and the exact equation (Boltzmann equation) governing 
transport phenomena can easily be derived, the solution to this equation 
is usually quite complicated. [1] 

Only in recent years have methods been developed to obtain exact 
solutions to the transport equation. These methods, while yielding exact 
results, are highly complex and require considerable computer effort. [2] 
For these reasons many problems in reactor physics and shielding design 
have been approached using the diffusion approximation to the transport 
equation. However, diffusion theory is expected to yield accurate results 
only when the assumptions of Pick's law apply to the problem. [3] 

Quite recently A. Leonard and 6. Garrettson [4] developed techniques 
to solve the neutron transport equation exactly for the class of problems 
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involving multidimensional neutron sources in a one-dimensional medium. 
Because some numerical data was available for a certain subclass of 
these problems, it was felt worth-while to develop a solution to the 
same subclass of problems using the diffusion approximation and to 
compare the results with the transport solutions. 

In this thesis, the Green's function for neutron diffusion in a 
homogeneous, isotropically-scattering slab, surrounded by a non reflecting 
medium, was determined using the diffusion approximation to the transport 
equation. In Chapter II, an analytical expression was obtained for the 
neutron density from a point source arbitrarily located in the slab. This 
was integrated to yield the neutron density from a pencil beam normally 
incident to the slab. A computer program was written to evaluate these 
expressions, and some of its features are described in Chapter III. In 
Chapter IV, the numerical results from transport theory are compared 
with those of diffusion theory. 

This thesis has two objectives. First, it was hoped that by comparing 
the numerical results of diffusion theory with those of transport theory, 
accurate estimates could be made of the region in which the diffusion 
approximation can be expected to yield satisfactory results. Second, it 
was hoped a systematic method could be found to choose the parameters 
used in the diffusion equation to yield optimal agreement between the 
results of diffusion theory and transport theory. 
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II. THE DIFFUSION THEORY SOLUTION 



A. THE EQUATION OF CONTINUITY 

Under the assumptions of one-speed, steady-state, and isotropic 
scattering, the equation of continuity in the diffusion approximation 
is 

Dv2$(_r) - (Eg - vE^)'i>(r) + S(_r) = 0 2.1 

where D is the diffusion coefficient; ${£) is the one speed neutron flux 
as a function of position, is the one-group absorption cross section, 

a 

E^ is the one-group fission cross section, v is the average number of 
neutrons produced per fission, and S(r^) represents an arbitrary source 
distribution function (those neutrons which have not yet suffered a 
collision). S(£) will later be considered a point source. Under the 
assumptions of Fick's law, the diffusion coefficient is given by 
D = e^/3e 2, where is the one group scattering cross section and 
^ is the total cross section [3]. 

The physical interpretation of eqn. (2.1) is as follows: 

-Dv2$(£)d2r represents the leakage rate out of a volume element, d^r, 
via streaming; E^$(r)d^r is the absorption rate in d3r;vE,$(r)d2r is 
the neutron production rate in d^r from fission; and S(r^)d^r is the 
production rate in d^r from any other source. 

To obtain the diffusion approximation to the transport equation one 
must make the following assumptions (Fick's Law): 

(1) The medium is infinite 

(2) The medium is uniform 

(3) There are no neutron sources in the medium 

(4) Scattering is isotropic in the laboratory coordinate system 

(5) The neutron flux is a slowly varying function of position . 

(6) The neutron flux is not a function of time. 
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The medium we consider is a homogeneous slab of thickness t < <» 
with isotropic scattering, and we solve the problem under the assumptions 
of one speed and steady-state. Thus assumptions (2), (4), and (6) are 
satisfied but (1), (3), and (5) are not. 

In order to obtain the diffusion equation (2.1), it was necessary to 
assume the diffusing medium to be infinite in all directions [3]. Hov/- 
ever, neutrons arriving from more than a few m.f.p. (mean free paths) 
contribute little to the neutron current density at a given point, so 
the diffusion equation is expected to be reasonably valid except within 
a few m.f.p. of the boundaries. 

The no source assumption, implying that all neutrons contributing 
to the neutron density are the result of collisions, is certainly not 
valid in the cases studied. The slabs we considered contained highly 
singular sources: the point source and the normal pencil beam source. 

However, since few neutrons originating from either source can be expect- 
ed to survive more than a few m.f.p. without suffering a randomizing 
collision, diffusion theory should produce satisfactory results at 
distances greater than a few m.f.p. from a source. 

The assumption of slowly varying flux is not satisfied near the 
sources or boundaries. In the derivation of Pick's Law, only first 
order terms in the Taylor expansion were retained [3]. Although second 
order terms will not contribute to the current density, terms containing 
third and higher order derivatives will make a contribution. Thus, it 
is necessary to restrict Pick's Law to regions where the second derivative 
of the flux does not change rapidly with distance. This specifically 
excludes regions near the sources we consider since both the normal 
pencil beam source and the point source (both represented using Dirac 
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delta functions) are highly singular. It should also be noted that 
since flux tends to vary rapidly in strongly absorbing media, it is 
also necessary to restrict Pick's Law to systems in which << E_. 

Such a restriction yields an average number of secondaries per collision, 



close to one [1]. 

From the preceding discussion, we can see that the diffusion theory 
approximations should produce satisfactory results in regions of the 
slab not too near boundaries or sources. Hopefully, our numerical results 
will indicate more precisely the region in which diffusion theory is 
accurate and the degree of inaccuracy outside this region. 

B. THE BOUNDARY CONDITIONS FOR THE STEADY STATE DIFFUSION EQUATION 

As is normally done in reactor boundary problems, the boundary 
conditions at the surface were chosen such that 

? Bn ■ ■ I 2-2 

where ^ is the normal derivative of the flux and i is the extrapolation 
length [3]. This is equivalent to saying that the flux vanishes at a 
distance i from the boundary. The value of i was chosen such that the 
solution to the diffusion equation (2.1) matched as closely as possible 
the more rigorous solution to the transport equation within the slab 
(away from the boundaries where the diffusion equation is valid). 

For a planar free surface, transport theory shows that i = 0.71 gives 
the best match, where is the transport m.f.p. [3]. This choice of 
£ should provide good agreement between transport theory and diffusion 
theory in the interior of the slab. However, since the diffusion equation 
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is not valid near boundaries, a solution obtained by this device will 
not give the correct density near the boundaries (e.g., the flux does 
not really vanish at a distance i outside the surface). 

An additional boundary condition is obtained in a subcritical medium 
(c < 1). The neutron flux decreases toward zero as the distance from 
the source increases. 



C. SOLUTION OF THE DIFFUSION EQUATION 

Consider the homogeneous slab of thickness, t, which has a total 
cross section, E, and which emits c secondaries per collision. The slab, 
infinite in both transverse directions, is surrounded by a vacuum or 
pure absorber. [See Fig. 1.] Under the assumptions of one-speed, steady 
state, and isotropic scattering, the equation of continuity in the dif- 
fusion approximation is given by (2.1). 

Since 4>(_r), the one-speed neutron flux, is a function only of the 
neutron density, n(r_), and the neutron velocity, v (assumed constant). 



$(r) = vn(r), 
eqn. (2.1) can be written 



Dv2 n(r) - (e^ - vE^)n(r) + = o. 



2.3 



Where xe(ji,t + a), and y, z ^(-“,“>). 

It is convenient here to introduce cylindrical coordinates, 

(x,£), £= (r,4>), and r = sj + z^ , which are illustrated in Fig. 2. 
Since it is assumed that the external source 



S(x,r.) 0, for Xz{z,t + £), (})e(0,2Tr) 

Y^-^oo 

the boundary conditions for our problem are 
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FIGURE 2. CYLINDRICAL COORDINATES 
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2.4 



n{0,^) = n(T + 2«.,?) = 0 for r e(0,°°), <{ie (0,2tt) 
and 

n(x,^) 0 for X e( £,t + «,)> <j> e(0,2Tr). 2.5 

/V 

r-H» 

To obtain the neutron density n(x,£) from any uncollided source, 
S(x,^), it is convenient to have the Green's function, G(x,f; x',7‘), 
written here as a function of cylindrical coordinates [ref. Fig. 2]. 

Then the formal solution to (2.3) for any source S{r)/\/ is given by 

7.1T ir- 

d(|)' \ dr' G(x,n x',r') ^ (x',£'). 2.6 

0 0 0 

The Green's function has the following properties: 

(1) G satisfies the homogeneous differential equation 

Dv2g - (z^ - vZ.)G = 0 2.7 

a T 

for all 

x e(t,T + «,), y, z e(-“,«>), and X 5^ x' , £ . 

(2) G satisfies the boundary conditions 

G(ji,n x',£') = G(t + x',f') = 0 2.8 

for all xe(x,,T + z) and y,z,y',z' [5] and G(x,‘^' x',r')^-> 0. 

(3) G satisfies the proper jump condition at _r = r ' . In cylindrical 
coordinates , [9] 

G(x,r; x' log |_r - I . 2.9 

r ^ r ' 




U 



Equivalently, G satisfies the diffusion equation 



Dv2G(r_',£') - - vi^)G(r^;_r') + 63(r_- r_‘) =0. 2.10 

with a point source, S(j 2 ) = 6^(r. - r'), where 



r' = (x',y',z'), x'e(t,T + t), and y', z'e 
6^ (r.-r‘) is the three-dimensional Dirac delta function given by 
63(j2 - jr ■ ) = 6(x - X') 6(y -y') 6(z - z') in cartesian coordinates and 

63(£ - 1 ^') = 62(F - £') 6(x - x') in cylindrical coordinates, where 

«2(r - ?■) = . V). [4] 

To solve for the Green's function one can find the eigenfunctions to 
the homogeneous form of equation (2.1) by separation of variables. Then 
6(r.; r') can be expanded in terms of these eigenfunctions, and the 
expansion coefficients can be chosen so that the jump condition is 
satisfied [7]. 

In cylindrical coordinates equation (2.10) is 



_1 + 

r ^ 91r 






9X 



j G(x,f; x'.r,') 






G(x,r;x' ,r') 



“ ■ ^ '^^(f - £') 6(x - x'). 2.11 

Assuming a homogeneous medium and a point source located at x‘ on the x 
axis, symmetry dictates that G(x,n x',0) is a function of r but not of 
4 >. In general then, G is a function of l£ - £' | rather than £ and 
separately. Defining 

g(x,lf - r' I ; x') = G(x,?; x' ,r') 

we obtain 
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r 



W 



_ 8 _ 

3r 



ax' 



•j g(x,r; x‘) - Q gCx.r'jX*) = - 6(x - x‘) 2.12 



where we have defined 



the coefficient 



Q = 



vE . 



2.13 



g(x-,f;x‘) satisfies the boundary conditions (2.8). For fixed x'e(0,t) 
and r'(0,«), the homogeneous form of equation (2.11), satisfied by the 
eigenfunction il»(x,r), is 



(If ' 

To find the eigenfunction tl»(x,r) we separate variables 
'|i(x,r) = X(x)R(f), 

and substitute into equation (2.14). This leads to the following 
separated eigenfunction problems: 



0 . FX = 0 



with boundary conditions 



and 



X(0) = X(t + 2a) - 0, . 



0 

dr^ r dr 



with the boundary condition 



2.15 



2.16 



* Note: For simplicity, the solution is forced to o at x = 0 and 
X = T + 2i,, rather than at x = - £ and x = x + £. For numerical work 
and comparison with transport theory, a change of variables, x = x - £, 
is introduced where x e(0,x) is inside the slab. 
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R(r) 0 

r->oo 



where 

a2 = Q + L2. 

Equation (2.15) has solutions of the form 

X„(x) = C, sin(L^x) t 4 cos (L^x) 

where the boundary conditions require that C 2 = 0 and that the eigen- 
values are 



By changing variables such that p - a^^r , where 'J ^ ^ 
equation (2.16) can be reduced to the modified Bessel s equation 



p, ^ . P R(P) = 0. 



2.17 



Since the solutions to equation (2.17) are the zeroth order modified 
Bessel functions I_(p) sod K (p), [6] the general solution is 



The boundary condition requires that C- = 0 since I (p) » . 

p-)-co 

Expanding the Green's function in terms of the eigenfunctions 



2.18 



we obtain 



g(x.|r-r'li x') = J, - f,' !)• 2-19 



Since 



X„,(x) R„(r) = ♦„{x.r) 
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satisfies the homogeneous equation (2.14) and the boundary conditions 
(2.8), the expansion (2.19) for G(x,n x',_r‘) certainly satisfies the 
conditions (2.7) and (2.8) required of the Green's function. It 
remains only to be shown that (2.19) satisfies the jump condition (2.9). 
This is obtained immediately since [6] 

^o^\\L - I) ^ t - \l- L'\ 

so that R^(i") satisfies the following differential equation [8]: 



3r^ F ""cTr 



m m 



R.(r) 



2ti 5 (f ) 
f 



2.20 



Thus equation (2.19) is the proper form for the solution. 

Substitution (2.19) into equation (2.10), and using the Fourier sine 
representation [8], 



cO 



6(X 






SI 



V. / niTTX 

"(ft 



2l 



2.21 



\ 



as well as equation (2.20), yields 



D(t + 2z) ( T + 2ii)’ 



where we have utilized the orthogonality property 

xni 

‘ mirx 



sin 



T + 2t 



s1n 



2£ 



dx = 



0 for m n 



0 



T + 2i 



for m = n. 



Therefore, recalling that G(x,n x',F‘) = g(x, jr. - I.' | > 
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ca 



g(x>lL-r' I ;x') = ^ ttD(t+2j^J (t+2£ ) ( t+2£ ) - 



fn- I 



2.22 



is the required Green's function for the point source problem. 

It would, of course, be possible to develop the Green's function 
for the normal beam problems in a like manner. However, this is not 
necessary since equations (2.22) and (2.6) represent the formal solution 



to (2.3) for any source. In particular, 
p(x,f;x') = 







g(x,|?-7'|; x‘) 5^(F)e dx' f'dr'd(ji' 

0 0 L 2.23 

is the neutron density in a slab from a pencil beam normally incident to 
the slab at x = y = z = 0, and this function satisfies equation (2.3) 

-zx' 



with 

V — 



. Thus p (x, |r-r' I ,x' ) is the Green's function 
for the class of problems involving neutron beams normally incident to 
a slab. Substitution of (2.22) into (2.23) and subsequent integration 
yields 

03 

d(x r-x'l = m ) l-(-l )'^exp[-z(T+2£)] ( . / mirx \ ^ 2 24 

PVX,r,x ; ^ ( (mir)^ + e( t+2a) ( x+2ii,y ^o^ m*^^' 

To facilitate comparison of our results with those of transport 
theory [1], we measure length in dimensionless units of mean free path 
(1 m.f.p. = 1/e) and shift coordinates in x and x' , x ->■ x - z, so that 
X £(0,t) lies inside the slab. Then the solutions (2.22) and (2.24) 
take the form 

CO 

g(x.nx') = C y K„tv> 2-25 
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and 



p(x,r;x') = C 







(imr)^ + (t+22,J 




lu-^l 

where all lengths are expressed in units of mean free path. For example, 



where Er = r/(l/E) is the radial distance expressed in mean free paths. 

For equations (2.26)-(2.27) , and from this point on, all distances are 
understood to be expressed in units of mean free path (e.g., e(t+2£) 

T+2i m.f.p.). The constants C and C will be adjusted to normalize the 
diffusion theory results to transport theory results at some chosen point. 




2.27 
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III. NUMERICAL TECHNIQUES 



A Fortran IV program (appendix 2) was written to compute numerical 
results using equations (2.25) and (2.26) for comparison with numerical 
data obtained by transport theory. Since both equations involve an 
infinite series, it was necessary to truncate the series at some point 
in order to obtain a solution. The purpose of this chapter is to 
justify the criteria used to terminate the series and also to explain 
the choice of the particular values used for the parameter Q/e2. 

A. TRUNCATION OF THE SERIES 

Both equation (2.25) and (2.26) are infinite series solutions. In 
order to obtain numerical output, it was necessary to truncate these 
series at some point in the summation process. That is, 



where we choose N large enough so that is less than some chosen 

e, where _ 



is the remainder term. The maximum absolute value of the AX (x) terms 

mm' * 



in both equation (2.25) and (2.26) is unity. Thus the Bessel function 
term dominates the series for a r >> 1. 




3.1 




3.2 



m 




3.3 



u 

factoring out noting that [6] 




3.4 
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with 



^ Tfzr 



for m >> 1 , 



yields 



00 




^ ) y k + N ( ® 



-fir_v 



3.4 



K=0 



This power series can be sunmed to yield a convenient estimate for an 
upper bound of the remainder: 



"Ni 



\ 

1 - e ^ t+2£ ' 



3.5 



In the computer program the series is truncated at N terms with N 
large enough so that 



K„(V> 

F E } 



3.6 



where we chose e sufficiently small to give the desired accuracy. Then 
(3.5) yields the upper bound 



^ < ^ 3.7 

'n 1 - exp[-f3,] 

for the remaining fraction. This upper bound is of the order of e except 

^ T "^’2 

when r << ~ — , so that the choice e = 0.0005 should give at least 
three place accuracy except for small r^. However, (3.7) represents a 
very conservative estimate of the upper bound, so we would expect good 
accuracy with this choice of e even when r is small. In fact, for "r = 0.05 
and N large enough so that 






0.0005, 



22 



numerical results indicate that ^ < 0.001, while (3.7) yields the 



R, 



N 



N 



^ < 0 . 01 . 



Conservative upper bound 

Therefore, on the basis of this discussion, it would appear that 
(3.6), with e = 0.0005, is a sufficient truncation requirement to insure 
three-place accuracy. 



B. . CHOICE OF THE PARAMETER Q/E^ 

In equations (2.25) - (2.27) the parameter Q/e^ appears, where Q is 
defined in eqn. (2.13). This parameter is a function of the one group 
material cross section of the slab. However, the one-group cross sections 
depend specifically upon the energy-dependent flux, which is not avail- 
able. Therefore, one would like to choose Q/e^ so that diffusion theory 
will yield results that agree as closely as possible with the results of 
one-speed transport theory. In this thesis it is demonstrated how Q/e^ 
might be chosen to accomplish this. 

For comparison, this parameter was chosen in two ways. The first 
choice assumes diffusion theory for a medium in which there is no fission. 
For the second choice, it is noted that as r increases, the transport 
theory solutions and the diffusion theory solutions decay like 
and I^q(“p^^)> respectively, where are the transport theory radial 
buckling modes [2]. To force diffusion theory to behave like transport 
theory for large r, the first buckling modes in diffusion theory and 
transport theory are chosen equal , = B-j . 

To obtain the expression for Q/e^ in the first case, recall the 
parameters 



and 
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c = 



E + vE. 
S f 



3.9 



defined in Chapter II. In the diffusion approximation [3] D 
and if there is no fission (e^ = D). Q/e^ reduces to 

Q = ^ = 3(1-c) 






3.10 



since c = + i.) in this case. 

Recall that the diffusion approximation demands a highly scattering 
medium, c c?1.0, so c = 0.9 is chosen for the numerical work presented 
here. For this value of c, eqn. (3.10) yields the value Q/e^ = 0.333. 

For the second case the parameter is obtained by comparing the 
transport theory and diffusion theory solutions. For large transverse 
distances from the source, |£-£' | >> 1 , a relatively simple transport 
theory solution is obtained [2] since the finite (M < «) point spectrum 
modes dominate the continuous spectrum modes: 






"s: 



n 



r (x) 
m 



r^(x') 

m 



0 [ e-li-’E- I] 



3.11 



fn- \ 

where m = 1,2,-**,M, are the transport eigenfunctions for the 

discrete eigenvalues 0 < B-| < 62 < ’ ” < £ 1 [2]. Using eqn. (3.4) 

it can be seen that the first mode dominates (3.11) for large radial 
arguments, |r-r'| >> 1: 



^ ^ ^1^^^ r^(x') + 0 [ 



G I 

“TTTTTTI-i' 

r-r 



Likewise, the first mode m = 1 dominates the diffusion theory solution 
(2.25) for !£-£' I >>1 . Therefore, if we choose the first radial buckling 
modes to be equal, a.j = 3-| , the two solutions will behave similarly for 

✓V/ ^ , 

|r - r' I » 1 . 
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With, this choice for (“r+^ 

expression is obtained for Q/r^: 



1/2 



, the following 



ef - (^+k) 



3.12 



Numerically this corresponds to Q/z^ = 0.2757, for c = 0.9 and t = 10 
and to Q/z^ = 0.1918, for c = 0.9 and t = 3. [2] 
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rv. DISCUSSION OF RESULTS AND CONCLUSIONS 



The following results were obtained using the computer program discus- 
sed in Chapter III. The cases presented are those for which the exact 
solution has been obtained by G. Garrettson [2] using transport theory. 

The purpose of this chapter is to compare the results obtained using dif- 
fusion theory with those obtained using transport theory and to evaluate 
the performance of diffusion theory for problems of this class. These 
comparisons were done for cases of homogeneous slabs of thickness t = 10 
and T = 3 m.f.p. (mean free paths) surrounded by a non-reflecting medium, 
under the assumption of one speed, steady state, and isotropic scatter- 
ing. All graphical results (Figs. 3 to 17) discussed in this chapter 
are found in Appendix 1. 

The first results presented. Figs. 3-6, are those for a pencil beam, 
normally incident at (0,0,0) to a slab of thickness x = 10 m.f.p., with 
a multiplication constant c = 0.9. The data has been normalized to the 
transport theory results at r = 4.4448, x = 5.1282 in order to facilitate 
the comparison. Figures 7-10 depict the results obtained from a pencil 
beam normally incident at (0,0,0) to a slab of thickness x = 3 m.f.p. 
with c = 0.9. Figures 7 and 8 were normalized at r = 2.01, x = 1.5, 
while Figs. 9 and 10 were normalized at r^ = 0.2, x = 1.5. Finally, Figs. 
11-16 represent a comparison of the point source solutions in a slab of 
thickness x = 10, with c = 0.9, for point sources located at depths of 
1 and then 5 m.f.p. 

In all but Figs. 11-13 the diffusion solutions are presented for both 
values of Q/z^ discussed in Chapter HI. In this chapter, and in Appendix 
1, Q/z^ is referred to simply as Q, since z = 1 in units of mean free path. 



26 



A. COMPARISON OF NORMAL BEAM SOLUTIONS 

It Is interesting to note the evolution of the neutron density from 
a pencil beam, N(x,r), vs, x e[0,T], as the radial distance, r, increases. 
For very small r (ref. Figs. 3, 9, and 10), the neutron density tends to 
have the same shape as the once-collided neutron density, which in turn 
resembles the uncollided density (a step function times an exponential). 
Transport theory predicts a very sharp discontinuity near the boundaries 
which is physically explained by the surface leakage [2]. Although dif- 
fusion theory correctly yields an exponentially decaying density within 
the slab,* it does not yield the steep gradient given by transport theory 
within about one half m.f.p. of either surface. The results of the two 
theories diverge partly because transport theory accounts for the 
preferential streaming toward the boundaries while diffusion theory does 
not. In fact, as previously discussed, the diffusion equation is simply 
not valid near boundaries. (Ref, Chapter II). 

For any x e(0,t), the diffusion theory solution appears to be most 
inaccurate for small radial distances, r < 1. This inaccuracy is apparent 
from both the normal beam and point source data (ref. Figs. 3, 9, 10, 11, 
and 14), and it arises partly because diffusion theory does not consider 
the preferential streaming from a singular source in the region of the 
source. Since Fick's Law is derived under the assumption that there are 
no sources in the medium (ref. Chapter II) we should not expect the 
diffusion results to be accurate near sources. In fact, within one m.f.p. 
of the source, the errors induced by using diffusion theory may be as 
large as 80 percent. However, beyond about three m.f.p. the two theories 
yield nearly identical results (if Q is computed using eqn. (3.12)), 



* Recall 
exponential . 



that a straight line on a semilog plot corresponds to an 
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B. COMPARISON OF POINT SOURCE SOLUTIONS 

For solutions to the point source problem, the most noticeable 
discrepancy betv/een the transport and diffusion data appears in the 
region within about one m.f.p. of the source (ref. Figs. 11 and 14). 

As previously mentioned, diffusion theory's inaccuracy near sources 
is partially due to its failure to consider the preferential streaming 
in -the region of the source. As the distance from the source is 
increased, the relative error decreases unless a boundary region is 
encountered. As expected, within one half m.f.p. of the boundary, the 
diffusion theory solution begins to diverge until a discrepancy of 
about 20 percent is noted at the edges. 

C. CHOICE OF THE PARAMETER Q 

It should be noted (ref. Figs. 5, 7, 8, and 16) that for large 
radial distances, r, the agreement between transport and diffusion 
theory can be greatly enhanced by the proper choice of the parameter 
Q (ref. Chapter III). The diffusion theory yields Q = ^ Tor a 
medium in which there is no fission. Hov/ever, it was expected that 




where B-| is the first transport eigenvalue [2], would yield better 

results. By examining both the normal beam and the point source data 

for the two choices of Q, the most optimum choice is apparently given 

by equation ( 3 . 12 ). This choice forces the diffusion solution to decay 

in r at the same rate as the transport solution (for large r^) because 

the first modes dominate for r>>l. The numerical data bears this out 

2 

since the Q = ~ ( t -(- ' 2y) solution is nearly identical to the transport 

solution for t ^ 3. Unfortunately, this choice of Q does not consistently 
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improve tKe results near ttie source (or near the boundaries) since the 
assumptions for pick's law are not valid in these regions and since the 
first mode does not dominate the solution for small r".. A plot of vs t 
is given in G. Garrettson's thesis for c = 0.9, and data for other c's 
is contained in an appendix [2]. 

A comparison of results in the x = 10 and x = 3 slabs (refs. Figs. 
3-10), indicates that diffusion theory is more accurate in thick slabs. 
This is to be expected since Pick's law is not valid near the boundaries, 
and in a thin slab the neutron is never very far from one of the two 
boundaries. 

D. CONCLUSIONS 

In general, one can say that the neutron density computed using 
diffusion theory can be expected to yield reasonably accurate results 
except as indicated in the shaded areas of Fig. 17. If one is not 
interested in the density within about three m.f.p. from the source and 
within about one-half m.f.p. from the edges, diffusion theory can success- 
fully be used. However, if one is interested in the density within these 
regions (e.g., the shaded areas in Fig. 17), it will be necessary to 
resort to the more rigorous and time-consuming techniques of transport 
theory to obtain accurate results. 

One shouod also consider that the diffusion theory solution requires 
very long to compute for small radial distances from the source, since 
the smaller ]£-^‘ | is, the more terms the truncated series includes. For 
example, computer runs which included r = 0.01 required as long as 18.5 
minutes to complete, which is comparable to the computation time required 
for the most difficult transport theory solutions. These runs, even 
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though lengthy, did not provide accurate data, and one must conclude 
that only the transport theory solution should be used to compute 
densities at very small radial distances. 

Lastly, the choice of the parameter Q is an important factor in 
obtaining agreement between the results of diffusion theory and transport 
theory in the regions where the diffusion equation is valid. The choice 
(3.12) significantly decreases the discrepancy between the tv/o results- 
at large radial distances from the source, e.g., f > 3 m.f.p. 

In shielding problems with beam sources, the effect of these 
discrepancies will probably be most noted when calculating the transmis- 
sion or reflection from the region of the beam. In these regions the 
diffusion results are most inaccurate since the emerging current from 
a region of the slab is a function of the neutron density in that region 
[3]. Diffusion theory should give satisfactory results for both trans- 
mission and reflection from regions not too near the beam. 

E. SUGGESTIONS FOR POSSIBLE EXTENSIONS 

For slabs sufficiently thin, no discrete poles exist for the transport 
equation [2]. In these cases, to choose a proper value for Q, an effec- 
tive mode of the transport equation would have to be found empirically. 
This could be done by making a semi log plot of the transport solution 
for fixed x, as a function of r (for large r), and using the slope to 
determine the effective decay constant. Mathematically, this would 
correspond to the smallest "effective" pole in the continuous spectrum 
of the transport operator. These effective poles should be investigated, 
not only empirically, as mentioned above, but also mathematically. 
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Further, it would be interesting to investigate the results obtained 

over a wide range of values of both t and c. One would expect diffusion 

theory to be less accurate for small t, because of leakage, and also for 
l + 

small c = ^ , which corresponds to a highly absorbing slab. The 

diffusion approximation is not expected to be valid unless << 

(Ref. Chapter II). 
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APPENDIX 1 - Numerical Results 



Appendix 1 contains graphical displays comparing numerical data 
obtained using transport theory with that obtained using diffusion 
theory. Figures 3-10 represent semi-log plots of the neutron density 
in a slab, as a function of normal depth, from a pencil beam normally 
incident to the slab, for various radial distances from the beam path. 
Figures 11-16 represent semi-log plots of the neutron density in a 
slab, as a function of normal. depth, from a point source within the 
slab, for various radial distances from the source. 

This data is discussed in Chapter IV. 
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APPENDIX 2 - Fortran IV Program Listing 



C STUDY OF NEUTRCN DIFFUSION 
C 

REAL+A NUB , NVAL, NUM, MQOB 
DIMENSION PdCOtlCO) , PT(100,1C0) 

DIMENSION XMESH (100) , R*^PSH (100) 

READ(5» 101 ) FPStXPRl ,TAU,C,D 
WRITE(6,101)EPS,XPRI,TAU,C,D 
READ(5»111)NVAL,PTVAL,P,Z,Y 
WRITF( 6, 11 1 )NVAL, PTVAL .R ,Z ,Y 
111 FORMAK 5F10o6) 

NR = 1 B 
NX = 35 

C THE 35 XMESH PTS ARE 

READ(5, ICCO ) (RMESH( K) ,K=1 ,NR) 

WRITE(6,100C) (RMESH(K) ,K=1,NR) 

READ (5,l?l) (X-^'ESH(J) ,J = 1,NX) 

•WRITP(6, 121 ) (XMESH( J) , J=1,NX) 

121 FORMAT <6Firo5) 

101 FORMAT ( 5F1C .5 ) 

ICOO FORMAT (5FlCo5) 

M = 0 

TEMP=0>0 

NUM = Ob C 

M0DB=0,0 

PI=3ol41593 

T=TAU+(2o^EPS) 

AR=0»0 
DCY=OoO 
ALPMR = Oo 0 
FUD=OoO 
CH=OoO 
A=OoO 
XSYN=0^0 
BESL = OoO 
TO=OoO 
NUB = 0.,0 
CK=OoO 
C2=OoO 
Z=7+EPS 
Y=Y + EPS 
XP=0«0 
XP=XPRI +EPS 
C 

C COMPUTE THE NORMALIZING CONSTANTS 
5 M, = M+1 
AB = M 

DKY=AB*( lo-(( (-lo )**M)*EXP(-T )))/(( AB '-^P I ) *'^ 2 + ( T*’f^2 ) ) 

C2=AB*PI/T 

A = R*S0RT(C + (C2*-"-2) ) 

XSYN=SIN(C2*Z ) 

CALL BESK( A,0, BK, I FR) 

IF ( lERo EQo C ) GO TO 200 
WRITE(6,201) lER 

201 format ( *0» , THE BESSEL ERROR I S ' t I 3 , ’ 1^**^'*' • ) 

200 CONTINUE 
BESL=BK 

TO=TO+( DKY*;‘XSYN’^BESL ) 

NUB = NUB +( 2o/T)*SIN(C2*Y)«=XSYN^BESL 
C COMPUTE SIZE OF BESSEL FN 
CK=RESL/TO 

1F(ABS(CK) -Co 0005) 125,5,5 
125 CONTINUE 

CCNN=NVAL/TO 

CCNP=PTVAL/NUB 

M,=0 

DO 3 K=1,NR 
AR=RMESH(K ) 

DO 4 J=1,NX 
X=XM.ESH< J ) 

AX=X+EPS 
120 M=M+1 
AM: = M 
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m 

ti 



non o on on non no oooo 



CCMPUTE DECAY TERM, DCY 

0CY = AM1' ( lo-( ( (-lo ) T + M) KEXP (-T )))/(( AM«'P + ( TY*2 ) ) 

CGMPUTW FRACTICN M*PI/T=C1 
Cl=AM*PI /T 



COMPUTE ALPMR 

ALPMR = AR*SORT( { D ) + ( C ) ) 

SIN X TER^ 

FUO= SIN(C1*AX) 

BESSEL FN term 

CALL Rf=SK( ALPMR,0, BK, I ER» 

IF { IERo EO., r ) GO To 202 
WRITE(G,201) lER 
202 CONTINUE 
M00B=BK 

TEMP=TEMP + (0CY*M00B*=FUD) 

CCMPLTE VALUE FROM POINT SOURCE 

NUM=NUM+( ( 2o / T) «:SIN(C1'^XP) *FUD*M0DB) 

C CHECK VALUE OF BESSEL FN 
C 

CH=M00R/TEMP 

IF ( ABS(CH)-0o0005 ) 14C, 120, 120 
lAO CONTINUE 

B( J,K)=CQNN*temP 
PT( J,K)=CONP"’=NUM 
M = 0 

NUM=0o0 
TEMP=0,00 
4 CONTINUE 
3 CONTINUE 

500 FORMAT! 'O' , 'TAU = •,F9,6,' ,C = 'tF9o6t' ,AND D = 
WRITF(6,500) TAU,C,D 

501 FORMAT(//'C' , 'THE* , 13, • X MESH POINTS ARE'/) 

502 FORMAT! 'O' , 10F12o6 ) 

WRITE!6,501) NX 

WRITF ! 6,502 ) !XMESH!J) ,J = 1 ,NX) 

503 FORMAT ! • 1 • ) 

WRITE! 6, 503 ) 

504 FORMAT !*o',«pt SOURCE AT DEPT H • , F 6/ • 0 ' , ' DENSITIES 
WRITE! 6, 504 ) XPRl 

505 FORMAT! 'O' , 'R = * , F 8o 4 , 4X , 1 P lOE 1 1 o 3 / ! • 0 * , 1 5 X , 1 P 10 E 1 lo 3 ) 
DO 105 K=1,NR 

WRITE! 6,505) R ME SH ! K I , ! PT ! J , K ) , J=1 , NX ) 

WRITE! 6,507 ) 

508 FORMAT !30!'»'I,'R =', F8o4,' FOL L OWS • , 30 ! ' t • ) ) 

105 CONTINUE 
WRITE! 6, 503 ) 

506 FORMAT !• O' ,' THE NORMAL BEAM DENSITIES ARE'//) 

WRITE ! 6,506 ) 

DO 106 K=1,NR 

WRITE!6,505) RMESH!K) , ! B ! J , K ) , J= 1 , NX ) 

WRITE !6,507) 

507 FORMAT ( '0 ' ) 

WRITE!7,508) RMFSH !K) 

509 FORMAT! IP 8 EIO 93 ) 

WRITE !7,509) !B!J,K), J=1,NX) 

106 CCNT INUE 
STOP 
END 



49 



I'ii 



BIBLIOGRAPHY 



1. Weinberg, A. M., and Wigner, E. P. , The Physical Theory of Neutron 
Chain Reactors , University of Chicago Press, 1958. 

2. Garrettson, G. A., Green's Functions for Multidimensional Neutron 
Transport in a Slab , Ph. D. Thesis, Stanford University, Stanford, 
California, 1965. 

3. . Lamarsh, J. R. , Introduction to Nuclear Reactor Theory , Addison- 

Wesley, Reading, Mass., 1966. 

4. Garrettson, G. and Leonard, A., "Green's Functions for Multi- 
dimensional Neutron Transport in a Slab," Journal of Mathematical 
Physics , V. 2, p. 725-740', February 1970. 

5. Friedman, B., Principles and Techniques of Applied Mathematics, 
Wiley, 1965. 

6. Handbook of Mathematical Functions; with Formulas, Graphs, and 
Mathematical Tables , National Bureau of Standards, 1964. 

7. Kraut, E. A., Fundamentals of Mathematical Physics, McGraw-Hill, 
1967. 

8. Carrier, G. F. , Krook, M. , and Pearson, C. E. , Functions of a 
Complex Variable , McGraw-Hill, 1966. 

9. Courant, R. and Hilbert, D. , Methods of Mathematical Physics , 

1st ed., V. 1, Interscience Publishers, 1965. 



50 



INITIAL DISTRIBUTION LIST 



No. Copies 



1. Defense Documentation Center 2 

Cameron Station 

Alexandria, Virginia 22314 

2. Library, Code 0212 2 

Naval Postgraduate School 

Monterey, California 93940 

3. Corps of Engineers 1 

Department of the Army 

Washington, D. C. 20315 

4. Asst. Professor G. A. Garrettson, Code 61 Gr 1 

Department of Physics 

Naval Postgraduate School 
Monterey, California 93940 

5. Major John H. Shimerda, USA 1 

USA CDC INS 

Ft. Bliss, Texas 79916 



51 



I 




ri 



UNCLASSIFIED 






DOCUMENT CONTROL DATA - R & D 

iSecurtty classification of title, body of abstract and indexing annotation must be entered when the overatt report is classified) 



1 OR>CiNATlNG A c T 1 V I T Y f Corpora /c aufrtor^ 

Naval Postgraduate School 
Monterey, California 93940 



2a. REPOR T SE CURI TV C L. A SSI F ( C A T i O 

Unclassified 



26. GROUP 



3 REPORT TITLE 



Neutron Transport from a Point Source in a Slab; A Comparison Between Diffusion 
Theory and Transport Theory 



4 DESCRIPTIVE NO TES fTypa o/ report inclusive da fesj 

Master's Thesis. September 1970 



5 AUTHORIS) (First name, middle initial, last name) 

John H. .Shimerda 



6- REPOR T D A TE 


7a. total no. OF PACES 


76. NO. OF REFS 


September 1970 


52 


9 


«a. CONTRACT OR GRANT NO 


9a. ORIGINATOR’S REPORT NUMBER(S) 


N/A 1 






6. PROJEC T NO 1 

N/A 


N/A 




c. 


96. OTHER REPORT NO(S) (Any other numbers that may be assigned 
this report) 




N/A 





10 DISTRIBUTION STATEMENT 



This document has been approved for public release and sale; its distribution 
is unlimited. 



n. SUPPLEMEN T A R Y NOTES 



N/A 



12. SPONSORING MILITARY ACTIVITY 

Naval Postgraduate School 
Monterey, California 



13. abstract 



The one-speed, steady state, diffusion equation was solved for a point source 
and for a normal pencil beam in a homogeneous, isotropically scattering slab. 
Numerical results obtained using diffusion theory were compared to available trans- 
port theory results for two slab thicknesses. 

This comparison demonstrated that the diffusion theory approximation to the 
transport equation will yield accurate results except within about one-half mean 
free path of a boundary and except within about three mean free paths of a source. 
The best agreement between diffusion theory and transport theory is obtained if 
the first radial buckling constant in the diffusion solution is chosen equal to 
the radial buckling computed using transport theory. 



DD. 



ro«. <473 

NOV 6S I “ / W 

S/N 01 01 -807*681 1 



(PAGE 1) 



53 



UNCLASSIFIED 



Security Classification 



1-31400 



UNCLASSIFIED 



Security Classification 



1 4 . 

KEY WORDS 


LINK A 


LINK B 


LINK C 


ROLE 


W T 


ROLE 


W T 


ROLE 


W T 


Neutron 
Radiation 
Diffusion Theory 
Transport Theory 
Boltzmann Equation 
Green's Functions 
Shielding 















DD ,?o?..1473 (BACK) 

S/N 0101 -807-6821 



54 



Security Classification 



A- 3 1 409 



2: 8 M A y 7 I 



1 8 9 6 i 



Thesis 

S473 

c.l 




Sh imerda 

Neutron transport 



from a point source 
in a slab; a compari- 
son between diffusion 
theory and transport 
theory. 

20 MAY 71 18961 



from a point source 
in a slab; a compari- 
son between diffusion 
theory and transport 
theory. 



Thesis 

5473 

c.l 



123340 



Shimerda 



Neutron transport 



fulfOli trj 




3 2 



isporl troni .1 poml source m 




768 001 95320 1 



OLIDLEV KNOX library 



